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A 3D DISCRETE MODEL OF THE DIAPHRAGM AND HUMAN TRUNK 

Emmanuel Promayon^ and Pierre Baconnier^'^ 

Abstract. In this paper, a 3D discrete model is presented to model the movements of the trunk during 
breathing. In this model, objects are represented by physical particles on their contours. A simple 
notion of force generated by a linear actuator allows the model to create forces on each particle by 
way of a geometrical attractor. Tissue elasticity and contractility are modeled by local shape memory 
and muscular fibers attractors. A specific dynamic MRI study was used to build a simple trunk model 
comprised of by three compartments: lungs, diaphragm and abdomen. This model was registered 
on the real geometry. Simulation results were compared qualitatively as well as quantitatively to the 
experimental data, in terms of volume and geometry. A good correlation was obtained between the 
model and the real data. Thanks to this model, pathology such as hemidiaphragm paralysis can also 
be simulated. 

Resume. Dans cet article nous presentons un modele discret 3D permettant de modeliser les mou- 
vements du tronc pendant la respiration. Les objets du modele sont representes par des particules 
physiques sur leurs contours. Une notion simple de force induite par des actuateurs lineaires permet 
de generer des forces au niveau des particules en utilisant un attracteur geometrique. Les proprietes 
elastiques et contractiles d'un tissu sont ainsi modelisees par des attracteurs de memoire de forme locale 
et de fibre musculaire. A partir d'une etude specifique en IRM dynamique, nous avons construit un 
modele de tronc simplifie comprenant trois compartiments : les poumons, le diaphragme et I'abdomen. 
Ce modele est recale sur la geometrie reelle. Nous confrontons les simulations obtenues aussi bien 
qualitativement que quantitativement, en terme de variation de volume et de geometrie. Une bonne 
correlation est obtenue entre le modele et les donnees reelles. Grace a ce modele nous montrons enfin 
que Ton pent simuler la paralysie hemidiaphragmatique. 



Introduction 

The diaphragm has two main roles: anatomicaUy it separates the thoracic compartment from the abdominal 
compartment and physiologically it is the main respiratory muscle. The action of this muscle is complex and 
depends mainly on its size, its shape, and its attachments and links to surrounding organs and skeleton. The 
human adult diaphragm is shaped like a dome: a central tendon originates the muscular fibers. Laterally the 
fibers are inserted on the 7th to the 12th ribs (see Fig. [1] left). During inspiration, the diaphragm contracts 
and the abdominal content plays the role of a lever resulting in an enlargement of the thoracic cavity. This 
enlargement generates a negative pressure inside the rib cage, drawing air into the lungs. When the diaphragm 
relaxes, the air is expelled, helped also by the elasticity of the lung and the tissues lining the thoracic cavity. The 
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abdominal compartment can be considered as incompressible during a given period of time (several minutes). 
Indeed, apart from a small gastric gas content, the abdominal cavity is filled with organs of quasi constant volume 
(blood volume variations are neglected) as all human tissues except lung. The stomach is commonly isolated 
from the remaining digestive tract by two closed sphincters, its gas content is then constant and considered 
incompressible in the range of observed gastric pressures. 

A model of the diaphragm and its surrounding structures can be used in two simulation fields: physiology 
and computer assisted surgery. It has to be geometric and kinematic as well as dynamic. If the simulated 
movements are produced by the model at a sufficiently fast rate, it can be used to predict the diaphragm and 
abdominal organ positions during respiration thus being able to drive an imaging device or a conformative 
radiotherapy. It is also important to be able to model specific diaphragm pathologies, such as hemidiaphragm 
paralysis, as they can highly alter the abdominal organ movements. 

Physiological studies of the respiratory system classically include volume and pressure variations. But as 
the diaphragm is not visible nor easily accessible from outside the body, studying the diaphragm deformation 
requires to use three dimensional medical images [21], either Computerized Tomography (CT) scan or Magnetic 
Resonance Imaging (MRI). Pettiaux and al [16] showed that CT scan allows satisfying 3D reconstructions of 
the diaphragm. Cluzel and al [5] and Craighero and al [7] shown similar results using MRI. There are few works 
dedicated to model the diaphragm muscle. Boriek and al. [3] used a Finite Element Method (FEM) membrane 
model to study the material behaviors but did not try to compare the model with experimental deformations. 
Kinetic modeling was also proposed in [12] and [15], using geometrical change to describe muscular actions. In 
respiratory physiology, the most famous model is a compartmental model where the rib cage and abdomen form 
two compartments and where an electric schema analogy is used to display the relationships between active 
and passive links. However it seems difficult to use this model to establish links between a given pathology 
and some local mechanical problems and to give 3D geometric information. Improvements of this model had 
recently been used to include planar geometry information [2]. Other models include computer graphics model, 
such as [23] and [24], and focus on computer graphics 3D animation rather than physiological realism. 
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Figure 2. An MRI image acquired using the Fast Field Echo-Echo Planar Imaging. The 
resolution is 256x256 pixels. Note the variation of contrast, and the artifacts due to the 
acquisition velocity. 

1. Method 

1.1. Available data 

Our current work is based on data acquired by a specific acquisition protocol [7]. A 1.5 Tesla MRI acquisition 
was performed using the Fast Field Echo-Echo Planar Imaging techniques. In conventional MRI each image is 
acquired in 5 seconds, has a resolution of 512x512 pixels and a thickness of 1 mm. In [7], the acquisition time 
was reduced to 227 ms in order to study the diaphragm deformations. The main drawbacks of this technique 
are the image resolution (256x256 pixels), the slice thickness (10mm) and the poor quality of the image due to 
movement and reconstruction artifacts (see Fig. [2] for an example of the acquired MRI data). Recently, using 
an enhanced protocol, MR images of three ventilated subjects were acquired. The respiratory volume and its 
variations were directly controlled by using artificial ventilation. From the MRI images, a post-synchronization 
process made possible the reconstruction of diaphragmatic surfaces (see Fig. [TJ right). 

1.2. Model 

As presented in the previous paragraph, the only available data are poor quality MRI images (see Fig. 
This yields to very strict specifications for the geometry and paramaters registration of the model. Another 
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requirement for the model, as stated in the introduction, is its ability to produce very fast and accurate 
simulation. Two different directions can be taken by researchers to model human soft tissues [9]: the classical 
biomcchanical approach and the computational discrete approach. The classical biomcchanical approach is 
mainly based on the Continuum Mechanics or uses its counterpart, the FEM. It offers the advantage of being 
based on a strong theoretical background. There are generally two kinds of drawbacks when one applies this 
method to computer aided medical or physiological simulation: the computation time cost and the difficulties 
to build complex assemblage where rigid, elastic and active structures are interacting. Computation time can 
be reduced in this context, even for material with non-linear constitutive law, by using recent derived methods 
such as [6,8,22] or alternative continuous models such as [1,10,13,20]. Finite Element Analysis is extremely 
interesting when one needs to understand and to know the consequences of a local deformation on the material 
stress and strain. However, in the present work, the aim is mainly to get an accurate, patient-specific geometry 
and dynamics. We need to know the consequences of the diaphragm contractions in terms of body structure 
displacements and deformations, i.e. at a higher scale than the many different materials composing the different 
organs and tissues. In the FEM, the extraction of physical parameters, such as Young modulus and Poisson ratio 
for linear constitutive law, is possible by measuring isolated tissue samples [11]. In vivo tissue characterization 
is essential because the mechanical behavior of soft tissues can differ significantly between in-vivo and ex- 
vivo conditions. Tissue characterization, as done in [14, 19], is nevertheless extremely difficult to perform on 
living tissue and/or in vivo, notably due to the tissue accessibility, the organ movements or deformations, and 
the need to sterilize the measurement devices. In this study the in-vivo measurement of the tissues rheology is 
impossible as no direct access to the organs is provided. Moreover the muscle activation function are not anyhow 
available. This means that a global optimisation process using the whole organ geometries and deformations 
during respiration has to be enough to fit the geometry and the physical parameters of the model. 

Our aims are to include multiple dynamic interactions and properties, to be able to produce real-time 
simulations, and to be able to fit the model only using the available MRI data. These aims justify the choice of 
a discrete approach. Previous works from the authors [17, 18], and more recently from Zordan and al. [23,24], 
used the same approach to build a visual simulation of the respiration. But in both cases, the simulation were 
not compared to patient-specific data nor even to experimental data. In this paper we propose to qualitatively 
and quantitatively compare the results of our discrete model simulation with the available MRI data. 

. To model living structures, we mainly need three different kinds of components: 

• rigid components (to model skeleton), 

• deformable components (to model soft tissues), 

• and active deformable components (to model muscles). 

In our model, these components are all derived from the same principle: a set of particles control the 
component surfaces, themselves organized using triangular facets. Each particle has a position, a mass and 
different properties depending on the kind of component it is part of. Accordingly, the particles in deformable 
components have an elastic property and the particles in active deformable components also have a contractile 
property. In essence, this is similar to a mass-spring network, but the elasticity is described using an original 
formulation, which allow better stability and control than mass-spring network, as shown in [18]. 

1.3. Dynamics 

Forces are exerted on the particles to generate displacements and deformations. Three kind of forces are 
needed : 

Force field: this kind of force is applied to all particles. At each time the force intensity and direction is 
known. This kind of force can vary depending of some mechanical or physical properties, e.g. the mass 
or the velocity of the particle. The gravitional force is an example of such a force. 

Focal force: this is a kind of force known in intensity and direction and applied at specific time of the 
simulation. For example it can be a force applied to a particle by one of its neighbor in a particular 
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type of interaction. This kind of force is also used to apply boundary conditions or to transmit user 
interaction in the model. 

Linear actuator force (LAF): this kind offeree depends on the internal state of the object, i.e. mainly 
on the geometry of neighboring particles. The intensity and direction of this kind of force is computed 
using local geometric or mechanical data and is generally updated at each time step. A LAF is used 
when a particle has to go toward an ideal position that minimizes a given function. Therefore, LAF are 
used to model elasticity and contracility. 
We introduce a LAF in order to minimize a given energy E. Whenever it is possible to define, for a given 
particle of 3D position P, a position PminE that is known to be minimizing E, a LAF can be used. Thus a 
LAF is simply a force that tends to minimize the distance |PPminE| ■ To express a LAF, we can use a simple 
expression such as: 

F = fciAF(PPmi„E) (1) 

where k^AF is a parameter of the particle, or of a whole components. LAFs can thus be seen as potential forces 
that tends to minimize a distance. LAFs can model any kind of forces that could be defined by a target position. 
PminE Can depend on geometry or on constraints. The most important and difficult part is to determine a 
correct expression for PminE, so that it approximates a local minimum of E. 

. Spring-mass network parameters are known to be difficult to find and adjust. Therefore, our model does not 
use a network of springs to link the particles. To model the elastic property of a particle we define a local 
elasticity memory [17]. The elastic property of a particle is simply its ability to come back to its original 
geometric configuration once deformed. To model this property each particle has a local coordinate system 
defined relatively to its neighboring particles. This local coordinate system is defined by three parameters: two 
angles a and (3 and one distance 7 [17], see Fig. [31 These three scalars are initialized at the rest configuration 
and are called ao, Po and 70. At any time, if a particle position verifies a — ag , /S = Po , and 7 = 70 , then 
the particle is at the rest configuration, thus the component is locally undeformed. Using this local coordinate 
system, we can compute at any time and for each particle a position using ao, /3o, 7o and the position of 
neighboring particles. This position is ensured to locally minimize the deformation energy. Using this position 
as PminE allows US to define a LAF to minimize this energy. 

An elastic component is defined by a contour where this particular LAF is applied to all the particles. The 
elasticity parameter is the stiffness keiasUdty used by the LAF. 

. Another LAF is used to model contractility. Once the contraction directions (muscular fibers) on the muscular 
component are defined, the position PminE of a particle at one extremity of a contractile components can be 
simply defined as being the particle position at the other extremity of the fiber. The LAF could then be 
activated by varying the kcontracuie — j'^it) -kmuscie coefficient during the simulation, kmuscie is constant. In 
order to activate the muscle contraction, A{t) mimics a muscle activation signal. A{t) can take all the values 
between [0..1]. When A{t) = 1 the activation is maximal, and when A{t) = 0, it is null. 

. To solve the system dynamics, at each time t, internal and external forces are computed, and the equation of 
motion are integrated, taking into account the local and global constraints. 

Note that a particle mesh can include any types of particle. For example, it is possible to have an elastic 
particle with muscular neighbors. This does not generate any computational problem. Each particle accumulates 
its internal forces (elastic forces or muscular forces) and corresponding reaction forces are then distributed to 
its neighbors, in order to verify Newton's second law, independently of their types. 

All the geometry and physical parameter are described using the PML language [4]. 

1.4. Constraints and loads 

Forces are not often sufficient to model complex behaviors. Constraints are added to maintain some conditions 
like non-penetrating area or incompressibility. Our algorithm considers constraints as non-quantified force 
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Figure 3. Local coordinate system use to define the shape function. For each particle three 
geometric parameters arc used; a, /3 and 7. Once determined at rest shape, they can be use 
to target the position that minimize the deformation. The considered particle is in red, its 
neighbors are in blue, G is the center of the neighbors masses, n the approximated normal to 
the surface around the particle. 

components: they are solved using a direct projection algorithm based on the gradient vector of the constraint 
function. 



Volume preservation. It is possible to handle the total incompressibility of a closed contour, and therefore 
to have a tighter link with real tissues. Improvements of the previously published method (see [17]) allows for 
real-time computation of this particular constraint and thus for any kind of triangulated surface. 

Volume preservation is an essential property of soft tissue modeling. The control of the volume is necessary 
in order to simulate both the incompressibility of some organs and to control the volume variation of other 
organs. 

Consider one object described by a triangular mesh of particles at the contour, in our model, the volume 
preservation constraint is applied to all these particles. Note that the triangular mesh can also be used for 
visualization. Let n be the number of particles of this triangular mesh. Let Pj denotes the positions of the n 
particles. Let Vb be the initial (rest shape) volume of the mesh and T^(Pi, • ■ • ,Pn) a function of the particle 
positions that gives the current volume value. If the volume-controled mesh is deformed during the simulation, 
our algorithm provides a fast and efficient way to preserve the inner volume while keeping the mesh shape 
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similar. Let P^, i G [1 • • • n], be the positions of the particles before the correction due to the volume constraint, 
that is to say just after the model forces have been summed and integrated on each particles. Our method is 
able to find the displacements to apply to each particle in order to correct the current volume. In order to find 
these displacements, the following system has to be solved: 

P. = P. + A|^,V»e [l-.-n] 
F(Pi,-- - ,P„) = Fo " 

where P^ are the corrected positions and A is the unknown scalar. is equivalent to a constrained 

corrective displacement that solve the volume-preservation constraint. Solving system ([2|) allows us to directly 
find a solution for the volume-preservation problem. By rearranging the equations, we can simplify system 
^ into an equation in A'^ which coefficients depend only on Pj and A^^. Compared to lagrangian methods 
(lagrangian multiplier and minimization algorithm) , our method exactly solves the constraint and is very fast as 
it is mainly the direct solution of a third degree equation. Note that this algorithm can also be used to control 
the volume variation of an object by modifying the targeted Vg during simulation. 



Boundary conditions. On top of all forces and the volume preservation constraint, any other boundary 
conditions can be applied, such as null or imposed displacement in any direction, and imposed forces. All this 
kind of boundary conditions are described using the LML language [4] , allowing for a dynamic change of the 
boundary conditions if needed during the simulation. 

1.5. Healthy and pathological diaphragm models 

Our discrete modeling framework was used to describe a simplified human trunk (see Fig. IH left). It includes 
three components: lungs, diaphragm and abdomen for a total of 113 particles. The lung is a passive area, and is 
only modeled geometrically to monitor the volume variation AV generated by the diaphragm contraction. The 
diaphragm is modeled using an elastic and contractile component. The abdomen is an elastic component. The 
model geometry was registered using an elastic matching algorithm to the geometry segmented and reconstructed 
from the MRI at the beginning of inspiration. 

The muscular fibers are defined on the model by selecting the particle that mimic the real muscular fibers: 
their direction is vertical and along the zone of apposition (see Fig. HI top right). The activation function is 
set to mimic the physiological signal (linear contraction for 2 seconds, important decrease for 0.5 seconds and 
then normal decrease for 2 seconds). To model the central tendon, we set keiasUdty of the top central part of 
the diaphragm as being twice as rigid as the other diaphragm areas (see Fig. lU bottom right). 

As the abdominal compartment can be considered as incompressible, the boundary conditions essentially 
consist in maintaining the incompressible constraint on the mesh defined by the diaphragm and the abdomen 
walls. A null displacement boundary condition is imposed to some particles to model the spine, the pelvis, and 
the top of the lung. 

The model was compared with the experimental data by studying the deformation during a respiratory cycle 
in different directions (Fig. [5]). We compared deformation in terms of the variations of, from the most significant 
to the least significant: pulmonary volume (AV), apposition zone length (AL), height of the diaphragmatic 
domes (AH), and transverse length (Ai?). 

We also simulated a pathological condition: an hcmidiaphragm paralysis. This was obtained by inactivating 
all fibers of the same side of the diaphragm. An additional modification was needed: the pathological model has 
to include the long term effect of the paralysis, namely the elasticity loss of the paralysed hcmidiaphragm. We 
thus set a different value of the elasticity parameter keiasUcUy for half of the diaphragm tissue. This pathological 
situation is known to lead to many ventilatory impairments among which a drastic decrease of inspiratory muscle 
efficiency, inducing a decrease of tidal volume (total volume displacement of each breath) and a paradoxical 
upward ("expiratory") movement of the paralysed hemidiaphragm during inspiration. 
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Figure 4. The healthy diaphragm model. The complete model using three compartments 
(lungs, diaphragm and abdomen) is shown on the left. The contraction fibers in the healthy 
diaphragm are represented as cylinders (top right). The central tendon of the diaphragm is 
shown on the bottom right. 
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% error 


-0.02 


-0.55 


-6 




6.57 









Table 1 . Qualitative comparisons between real data and model. 



1.6. Estimation of the model parameters 

In order to estimate the model parameters, the only available data were the pulmonary volume variation Ay. 
The main advantage of our model is its very fast computation time and its reduced number of parameters. In the 
trunk model, only two parameters are to be estimated: keiasUdty and kmuscie- An optimization algorithm based 
on an "analysis by synthesis" strategy was elaborated. It consisted in a four step loop: (1) assume a given set 
of parameters, (2) build and simulate a respiratory cycle using the model, (3) compare the provided simulations 
with the respiratory volume measurements in the least square sense, (4) from this comparison deduce better 
values of parameters in order to improve the simulation/measurement fit. This loop was continued until the 
comparison carried out in (3) gives satisfactory results. 
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Figure 5. Comparison measurements. Deformations are measured in the zone of apposition 
(AL), the height of the diaphragmatic domes {AH) and the transverse direction (Ai?). 

2. Results 

The simulation of a complete respiratory cycle only takes 1.50 seconds on a Pentium Xeon 5140 at 2.33 Ghz, 
i.e. a frame rate of approximately 3000fps. Comparisons between the model and the real data are presented in 
Table 1. Qualitative geometry comparisons were also made between the surface of the diaphragm in the model 
and the reconstructed diaphragmatic surface at the end of inspiration (see Fig. [6]). 

The model is able to simulate the hemidiaphragm paralysis pathology. Comparisons between healthy and 
pathological diaphragm can easily be observed in 3D (see Fig. [7]). 

We also can see significant differences between lung volume displacements (see Fig. [5]). 

3. Discussion 

3.1. Healthy diaphragm 

The model was able to reproduce an accurate volume variation AV and piston-like deformation AL. The 
model deformation compared to real deformation measured on a subject was as well qualitatively satisfying. 
Another very important point, especially when considering the application of the method in computer assisted 
medical intervention and physiological studies, is the fast computation time: the simulation is about four 
times faster than the respiratory cycle it is simulating. This result leaves some space for the improvement and 
enhancement of our model. 

On the other side the simulation are far from correct when AH and AR are compared. These differences are 
probably due to an over simplification in the discretization of the diaphragmatic zone and to the model itself, 
which does not take contact and friction into account. 

Although discrete and very simple, this model efficiently reproduced the complex movements of breathing. 
The major drawback of this model is that being discrete, it is not possible to compute the extract strain 
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Figure 6. Comparisons of the model diaphragm geometry (coarse mesh) and reconstructed 
diaphragmatic surface. The model was initially deformed to match the reconstructed shape at 
beginning of the respiratory cycle (end expiratory position) (left). After the simulation of a res- 
piratory cycle, the simulated deformation are superimposed with reconstructed diaphragmatic 
surface at the end of inspiration (right). 




Figure 7. Comparison between a healthy diaphragm model and hemidiaphragm paralysis 
model using the same activation function at end inspiration (during exercise breathing). 

and stress on the different components. As these values on the in vivo diaphragm are not obtainable by any 
technique, the choice of a continuous model does not seem to be crucial. 

3.2. Hemidiaphragm 

We observed all the clinical consequences of the simulated pathology. The paradoxical behaviour of the 
paralysed hemidiaphragm (upward displacement during inspiration) is evidenced on the 3D simulation (Fig. [71 
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Figure 8. Comparison of the volume variation during exercise breathing generated by the 
healthy diaphragm model (continuous line) and the hemidiaphragm paralysis model (dashed 
line). 

right). The difference between volume displacements in the healthy and pathological diaphragm amounts to 
what is typically measured (a 50% shortening, T. Similowski, personal communication) in clinical results (Fig. 
[5]). 

3.3. Future works 

Once the k^iasticity parameter is set for a given subject, the main advantage of this optimization technique 
is that we can directly and quickly adjust kmuscie depending on AV. This can lead to a real-time prediction of 
the diaphragm position during breathing, considering only one medical image at the beginning from which the 
diaphragm geometry can be registred. This could be used for example during conformative radiotherapy. Future 
works on the diaphragm model will include testing and validating other breathing situations and comparisons 
with other subject data. To overcome the differences noted for AH and Ai? we also plan to add the rib cage 
and its cartilage components (this work has just started, see Fig. [9]). 
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